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1.  INTRODUCTION 

The  compressive  behavior  of  coiiq)osite  materials  has  been  studied  extensively  during 
the  past  three  decades  and  a  review  of  literature  on  the  subject  is  beyond  the  sc(^  of  this 
paper.  Substantial  listings  of  references  on  the  subject  can  be  found  in  the  articles  by  Shuait 
(1985),  Camponeschi  (1991),  Guynn  et  al.  (1992)  and  Piggott  (1993).  Suffices  to  say  that 
the  compressive  response  of  composites  was  found  to  depend  on  the  prtq)eTties  and  response 
of  the  constituent  materials  and  on  the  fiber  volume  fraction.  As  may  be  expected, 
compressive  strength  is  sensitive  to  imperfections. 

The  essential  novel  feature  in  the  present  work  is  the  incorporation  of  random  fiber 
spacings,  as  commonly  encountered  in  composites,  into  a  model  for  their  compressive 
behavim*.  The  main  consequence  of  the  foregoing  feature  is  that  it  predicts  a  response  which 
involves  the  emergence  of  highly  concentrated  lateral  forces  on  the  fibers  simultaneously 
with  micro*buckling.  These  lateral  forces  are  a  most  likely  cause  for  the  development  of 
kinks.  One  of  the  outstanding  issues  regarding  the  compressive  response  of  composites  is 
that  the  common  methodology  for  predicting  con^ressive  failure  stems  hnom  considerations 
of  buckling  and  stability,  wtule  most  failed  specimens  exhibit  localized  kink  bands  which 
span  the  thicknesses  of  the  test  coupons.  It  seems  that  all  other  models  address  micro- 
buckling  and  micro-kinking  exclusively  of  each  other,  and  can  thus  be  grouped  accordingly: 

(1)  Models  which  consider  buckling.  These  include  the  work  of  Rosen  (1965), 

which  seems  to  be  the  first  article  on  compressive  failure  of  composites.  Considering  "shear 
-mode  buckling"  that  model  predicted  a  failure  stress  Oot  =  0^/(1  -  Vf),  where  Gm  is  the 
shear  modulus  of  the  matrix  and  Vf  the  fiber  volume  fraction.  That  prediction  is  inadequate 
for  two  reasons:  (a)  it  gives  which  is  several  times  higher  than  experimental  values,  (b) 
the  relation  Odt  ~1  /  (I  -  Vf)  contradicts  experimental  observations  which  show  that  ot:R 
grows  linearly  with  Vf  (at  least  up  to  Vf  =  0.55)  (c.g.  Piggott  and  Harris  (1980),  Morley 
(1987)). 

Several  modifications  to  Rosen's  model  were  introduced  subsequently.  Primarily, 
these  modifications  considered  non-linear  shear  response  of  the  matrix  and  initial  fiber 
waviness  (e.g.  Wang  (1978),  Lin  and  Zhang  (1992),  Guynn  et  al.  (1992),  Highsmith  et  al. 
(1992)  and  others  listed  in  the  aforementioned  review  articles).  Additional  modifications 
included  the  incorporation  of  fibers'  shear-deformation,  such  as  by  Davis  (1975),  or  the 
accounting  for  large  deformations  of  the  fibers  by  Yin  (1992).  Though  the  latter  model 
stems  from  a  buckling  formulation,  it  is  worth  noting  that  it  proposes  a  criterior  for  kink 
fOTination,  which  occurs  when  fibers'  curvature  attains  a  critical  value. 
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(2)  Models  which  consider  the  a-priori  existence  of  kinks:  These  include  works  by 
Evans  and  Adler  (1978),  Hahn  and  Williams  (1986),  and  Budiansky  and  Fleck  (1992). 

The  compressive  response  of  muld-diiecdonally  reinforced  laminates  such  as  Shuait 
(1989)  and  of  cylindrical  shells  such  as  Blake  and  Staibuck  (1993)  is  beyond  the  scope  of 
this  article.  Suffices  to  say  that  these  complex  circumstances  activate  various  modes  of 
failure  which  do  not  occur  in  the  uni-directional  case  considered  herein. 

In  all  the  above  works  the  riber  reinforced  composites  were  viewed  as  lamellar 
regions  which  consist  of  riber  and  matrix  layers  as  shown  in  Figure  la.  It  should  be  noted 
that  several  investigators  (Sadowsky  et  al.  (1967),  Hemnaim  et  al.  (1967),  Lanir  and  Fung 
(1972)  and  Greszczuk  (1975))  considered  fibers  of  cylindrical  geometry.  All  those  works 
assumed  linear  elastic  behavior  of  fiber  and  matrix  materials. 

In  addition  to  random  fiber  spacing  the  current  model  includes  initial  fiber  waviness 
and  considers  the  non-linear  shear  stress-strain  response  in  the  matrix.  The  fibers  are 
assumed  to  deform  in  accordance  with  classical  beam  theory. 


2.  FORMULATION  AND  RESULTS 

Consider  a  uni-axially  reinforced  composite  which,  following  Rosen  (1965),  is 
represented  by  a  two-dimensional  layered  array  as  shown  in  Figure  1(a).  Let  x  and  y  denote 
Cartesian  coordinates  in  directions  parallel  and  transverse  to  the  layers,  and  designate  by  2h 
the  thickness  of  a  "fiber  layer"  centered  within  a  composite  layer  of  thickness  2c.  Conse¬ 
quently  we  have  Vf  =  h/c  and  =  (c-h)/c,  where  Vf  and  Vm  are  fiber  and  matrix  volume 
fractions,  respectively. 

We  focus  attention  on  the  "shear  mode"  of  buckling  (Rosen  (1965),  Garg  et  al. 
(1973)),  where  all  fibers  buckle  in  phase.  Then,  following  Rosen's  premises  (1965)  for 
high  perfcnmance  composite  material  systems,  we  assume  that  the  external  compressive  load 
N  is  borne  entirely  by  the  fiber  region,  which  is  modelled  as  a  Bemoulli-Euler  beam,  while 
the  matrix  responds  in  shear  only.  Consequently,  we  have  the  following  familiar  expression 
for  YSy,  the  shear  strain  in  the  matrix: 

ym  a— I  -dvl  (1) 

1-Vf  dx  ^  ^ 

In  equation  (1)  v^  denotes  the  lateral  displacement  ( in  the  y-direction )  of  the  fiber.  In  view 
of  the  assumption  of  Bemoulli-Euler  theory,  v^,  and  thereby  also  y^V,  depends  only  on  x. 

In  addition,  we  consider  a  micro-buckling  length  L  and  and  initial  fiber  waviness 
v^{x)with  periodicity  of  2L.  In  the  sequel,  we  let  v  ^  (x)  =  5o  cos  (tcx/L),  though  this 
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specific  choice  is  not  essential  to  our  method.  Finally,  in  anticipation  of  the  ciicumstances 
which  emerge  due  to  non-uniform  fiber  spacings,  we  denote  by  q(x)  the  distributed  lateral 
load  on  the  fiber.  See  Figure  1(b). 


Considering  non-Unear  shear  resptmse  of  the  matrix,  we  write 

t!!V  =  G?F(Y5V)  (2) 

where  the  function  F(YSy)  expresses  that  non-linear  shear  behavior  of  the  matrix 
scaled  by  the  initial  shear  modulus 

The  longitudinal  strain  in  the  fiber,  e(,  under  the  combined  effects  of  compression 
and  bending  is  given  by 


In  equation  (3),  u^  denotes  the  fiber  displacement  in  the  direction  of  x. 
Consequently,  the  axial  displacement  at  x  =  L/2  is  given  by 


(3) 


dx+i^ 


(4) 


As  can  be  noted  fiom  equations  (3)  and  (4)  the  hypothesis  that  fibers  deform  in-phase 
in^lies  that  u^  and  v^  ate  commcxi  to  all  fibers  regardless  of  their  spacing.  On  the  other  hand, 
equations  (1)  and  (2)  state  that  the  support  provided  by  the  matrix  varies  with  the  fiber  vol¬ 
ume  firaction  Vf.  These  observations  imply  the  existence  of  lateral  loads,  q  =  q(x),  which 
enforce  a  common,  in-phase  deformation  of  all  fibers  in  the  case  of  non-uniform  spacing. 
To  en^hasize  their  dependence  of  the  spacing  c,  we  shall  write  q  =  q(x,c). 

Consider  an  individual  cell  of  width  2c.  The  principle  of  virtual  woiic  yields 


I  o{5ej[dVf  + 1 


-[  ..(X,. 


5ej[dV^+ I  q(x,c) 6vfdx +  N 5A  =0 


(5) 


Substitution  of  expressions  (l)-(4)  into  equatitm  (5)  and  employment  of  integrations  by  parts 
yield  the  following  field  equation  and  boundary  conditions  for  each  individual  cell: 
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with 

<t^  =  0,  <^4=0  at  x-0 

dM  I  <’> 

vt  =  0,  ^.0  at  x-L 

Note  that  in  view  of  the  non-linearity  of  F  in  its  argument  equation  (6)  is  a  non-linear 
differential  equation  for  v^. 

Turning  to  the  case  of  random  fiber  spacing,  let  p(c)  denote  the  probability  density  of 
the  cell  dimension  2c.  Obviously 


f 


p(c)  dc  =  1 

In  the  present  circumstance  the  principle  of  virtual  work  gives 

oC5eUV^  +  J  t5V5Yi!VdV”-|  q(x.c)  5vf  dx  +  N  5A 


(8) 


1  "'''(I 


dc=0  (9) 


Furthermore,  in  the  absence  of  external  lateral  loads,  equilibrium  in  the  direction  of  y 
requires 


(10) 


Integration-by-parts  of  equation  (9),  upon  expressing  all  variations  in  terms  of  5(^),  gives 

dx 

the  following  field  equation  for  v^ 
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The  boundaiy  conditions  for  the  case  of  randomly  spaced  fiber  remain  the  same  as  those 
given  in  equations  (7).* 

It  is  advantageous  to  further  reduce  the  order  of  the  differential  equation  given  in  (1 1) 
and  express  it  in  a  non-dimensional  form  in  terms  of  the  following  non-dimensional 
parameters: 

x  =  i.  Y.^.  e=^  (12) 

In  addition,  the  probability  distribution  function  p(c)  can  be  converted  to  a  probability 
oistribution  function  p  (Vf). 

In  view  of  expressions  (12),  the  non-dimensional  form  of  Equation  (1 1)  reads 


0-[  p(Vdo»(l-V0F(j^]dVr  +  X2Y  =  -X*Y„  (13) 

where 

(v2  - _ 21l-  ■  -  NL^  Y  =  -  EJi  sin  xX  (141 

Vf(l-Vf)  Elf  ’  Elf  ’  ®  eJtsinxA  (14) 

The  boundary  conditions  which  accompany  the  second  order  non-linear  differential  equation 
(13)  are 

In  the  case  of  uniform  fibe^  spacing  equation  (13)  reduces  to 

^-a2(l-Vt)F(^) +2.*Y  =  -2,2y„  (16) 


Note  that  for  the  linearly  elastic  case  with  uniformly  spaced  fibers  F(Y/(1  -  Vf))  =  Y/(1-Vf) 
and  equation  (16)  takes  th<*  simple  fnm 

^-a^Y  +  X^Y  =-X^Yo 
dX2 


*  In  view  of  equation  (10)  it  was  possible  to  derive  differential  equation  (11)  which  is  one  (xder  lower  than 
that  given  in  equation  (6).  The  lower  order  equation  (11)  enables  the  determination  of  the  lateral  displacement 
vf  to  within  a  rigid  translation,  which  is  of  no  relevance  to  the  failure  mechanisms  considered  in  this  work. 
An  additional  integration  of  expression  (11)  with  req)ect  to  x,  further  reduces  the  order  of  the  differential 
equation,  leading  to  a  solution  which  incorporates  an  indeterminate  rigid  body  rotation. 
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with  the  solution 


Y  =  — - sinjcX 

X2_K2_a2 

This  conesponds  to  the  buckling  load  predicted  by  Rosen  (1965),  namely  s  +  a^. 
Note  that  the  above  result  assumed  that  the  magnittide  of  die  lineariy  elastic  shear  strain  in  die 
matrix  is  not  limited  by  any  ultimate  or  plastic  level. 

However,  if  one  ctmsiders  a  linearly  elastic  matrix  response  followed  by  an  ideally 
plastic  defonnation  at  yjy  s  then  plastic  yield  begins  at  X  =  1/2  and  the  onset  of  plastic 
deformatitm  is  found  to  occur  at 

^  ent  (l-VdYp  ^ 

The  above  result  agrees  with  the  value  obtained  by  Steif  (1988)  for  the  slippage  initiation 
load,  beyond  which  the  matrix  no  longer  supports  the  deformed  Hbers. 

rasa  t  r  T  Tniformlv  Spaced  Fibers  with  Bi-linear  Shear  Modulus  of  the  Matrix 

Consider  a  bilinear  shear  stress-strain  response  of  the  matrix  material,  given  by  the 
following  expression  for  F(y|^) 


G?F(Y,y)  = 


^G?'(Yxy-Yy)+G?yy 

Ge*Yxy 

iG5‘(Yxy  +  Yy)-G?Yy 


if  Yxy>Yy 
if  -Yy<Yxy<Yy 
if  Yxy<-Yy 


(18) 


In  equation  (18)  is  the  strain  level  where  the  slope  of  the  bi-linear  stress-strain  diagram 
changes  from  an  initial  value  GS'  to  the  strain-haidetung  value  G^.  It  will  be  shown  in  this 
section  that  the  buckling  associated  with  the  response  expressed  in  equations  (18)  can  be 
handled  analytically. 

For  loads  that  correspond  to  X^  which  exceeds  X2(Yy)  in  equation  (17)  the  shear 
response  of  the  matrix  will  follow  the  H-linear  stress-strain  relation  over  a  region  4<X<l/2, 
but  will  still  remain  linearly  elastic  within  the  central  region  0<X<^.  Obviously  ^  decreases 
with  increasing  X^.  Substitution  of  expressions  (18)  into  equation  (16)  gives 
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(19) 


-aiY  +  X*Y»-X*Yo  at  0<X<^ 

dX* 

^ -ajY-»-X*Y=-X^Yo-p^  at  5<X<^ 

where  Oe  and  Op  defined  according  to  (14)  widi  shear  naoduli  GS*  and  G^,  respectively,  and 

.2-  2h  (GP-Gg)L» 

'*  Vt(i-Vt)  m, 

The  boundary  and  continuity  conditions  associated  with  equations  (19)  are 


-  0,  Y(0)  -  0,  .  Yft  ).  X  Y(y  —  V„Y,  (20) 

The  above  conditions  conespond,  respectively,  to  the  vanishing  of  the  moment  at  xsL/2,  and 
of  the  shear  at  xM),  the  ccmtinuity  of  shear  and  moment  at  x^L  and  the  requirement  that,  by 
hypothesis,  I  I  =  Yy  at  x^L.  The  tive  conditions  given  in  equations  (20)  determine  the 
four  unknowns  associated  with  the  two  second  cnider  differential  equations  (19),  as  well  as 
the  yet  unknown  location  ^ 

Note  that  the  solution  for  Y  determines  the  displacement  v^  to  within  arbitrary  rigid 
translations  and  rotations,  which  are  determined  from  the  requirement  of  continuity  of  v^  and 

at  x=£L,  as  well  as  vf(0)  =  0  and  ^  =  0  at  X5=L/2. 
dx  ^  dx 

The  solution  to  equations  (19)  reads: 


forO<X<^ 


Y  ^  sinh  tCeX 
Sinh  KeS 


I  - sin  7t^  +  (1-V^y  \ 

\  X"*  -  7t2  -  (xi  j 


+ - £EA - sintcX  (21) 


for^<X<l/2andX^<(x5 
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•»  ,  - sinicX- 


_JL> 


while  for  4<X<l/2  and  X^  > 


Y^) 


cosh  Kp(l-2X)/2 
cosh  i^l-24)/2 


,  ■  sin  it4  +  (1-Vf)yy  - 

X^-k^-(4 


(22b) 


+ _ eaX? _ 

X^-jc2-a^ 


sin^X- 


X^-aJ 


In  the  above  equations  Ke  -  Vet?  -  X^  and  Kp  =  Vlo^  -  X^l. 

Equadons  (21)  and  (22)  match  all  the  conditions  (20)  except  the  continuity  = 

dX 

The  latter  condition  yields  a  characteristic  equation,  upon  differentiation  of 

equations  (21)  and  (22),  which  relates  the  position  of  ^  to  the  load  parameter  X^.  This 
characteristic  equation  must  be  solved  numerically,  with  the  physically  meaningful  solution 
emresponding  to  the  lowest  value  of  X^. 

In  our  computations  we  utilized  the  constituent  properties  reported  by  Guynn  et  al. 
(1992)  for  AS4/PEEK  at  21*  C  Accordingly,  we  took  Ef  *  67  GPa,  L  »  330  pm  and  So  = 
1.65  pm  and  Vf  =  0.6.  For  purposes  of  coiiq>arison  we  also  considered  additional  values  of 
Vf  in  the  sequel  The  non-linear  shear  stress-strain  response  was  approximated  by  a  bi-linear 
leladonship  with  G?  =  1.3  GPa,  G^  =  0.33  GPa  and  7y  =  4.2%. 

The  resulting  stress-deflection  curves  ate  shown  in  Figure  2  for  various  values  of  Vf. 
The  symbols  on  those  curves  correspond  to  load  and  displacement  values  at  onset  of 
departure  from  linearity  in  the  shear  stress-strain  response  of  the  matrix.  Such  departure 
occurs  when  I  I  =  at  X  =  1/2.  Note  that  when  Vf  =  0.9  the  conqwsite  can  carry  com¬ 
pressive  loads  which  exceed  the  level  which  cause  departure  from  linear  matrix  response. 
However,  for  Vf  =  0.3  and  Vf  =  0.6  the  stress-deflection  curves  exhibit  the  so  called  "finite 
disturbance  buckling  behavior,"  resembling  the  buckling  of  cylindrical  shells  under  uniaxial 
conqnession  or  spherical  shells  under  external  pressure  (Simitses  (1976)).  It  is  interesting  to 
note  that  for  Vf  «  0.3  and  Vf  =  0.6  the  cusps  in  the  stress-deflection  curves,  which 
correspond  to  maximal  load  levels  prior  to  buckling,  occur  at  magnitudes  just  above  those 
which  cause  I  ysy  I  =  Yy  at  X=l/2.  It  is  obvious  that  the  dieoretically  predicted  cusps  for  Vf = 
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0.3  and  0.6  cannot  be  realized  experimentally.  Under  load  controlled  tests  the  maximal  loads 
will  be  followed  by  total  collapse  and  under  displacement  controlled  circumstances  the 
specimen  would  snap  through  to  the  lower  load  levels  along  the  vertical  dashed  lines  shown 
in  Figure  2. 

Further  insight  into  the  compressive  response  predicted  by  the  solution  to  equations 

_  ^  ^ 

(19)  and  (20)  is  provided  in  Figures  3  and  4.  The  dimensionless  length  ^  (^=1/2  -  4  )  of  the 
regions  where  the  matrix  shear  strain  I  i  exceeds  the  linear  elastic  limit  Yy  is  plotted  vs. 
the  applied  compressive  stress  Oc  in  Figure  3  for  fiber  volume  fractions  Vf »  0.3, 0.6  and 

0.9.  Note  that  Oc  increases  monotonically  with  ^  for  Vf  s  0.9,  but  decreases  (after  very 
slight  initial  amplifications )  for  Vf = 0.3  and  0.6. 

The  variadon  of  the  matrix  shear  strain  with  the  dimensionless  distance  X  along 
the  fiberAnatiix  interfaces  is  shown  in  Figure  4  for  Vf  =  0.6.  The  four  curves  in  that  figure 
correspond  to  disdnet  levels  of  non-dimensional  load  X.  The  top  curve,  with  X  =  23.10 

represents  typical  linear  elastic  results,  with  I  y  ^  I  <  yy  for  all  X  and  thereby  also  ^  =  0.  In 
this  case  we  obtain  a  sinusoidal  variadon  of  yiny  which  agrees  with  earlier  results  (Wang 
(1978),  Lin  and  Zhang  (1992)),  namely  yij\  =  A  sinnX  with  A  =  enX^  /  [(l-Vf)(X2.ji2.(x2)] 
The  foregoing  sinusoidal  variadon  persists  until  the  onset  of  inelastic  response  at  X=l/2 
which  occurs  at  X  =  Xy  =  30.79.  This  result  is  shown  by  the  dashed  line  in  Figure  4.  The 
maximal  value  of  the  coo:q)ressive  load,  associated  with  X  =  Xmax  =  30.81,  corresponds  to  an 
inelastic  zone  of  dimensionless  length  ^  =  0.03.  In  this  case  the  variation  of  yS'y  with  X, 

shown  by  the  dotted  line  in  Figure  4,  is  no  longer  sinusoidal.  Beyond  ^  =  O.OS  values  of  X 
decrease  while  d/L  increase  according  to  Figure  2.  A  typical  circumstance,  corresponding  to 
^  =  0.1  and  X  =  30.23,  is  shown  by  the  solid  line  in  Figure  4. 


As  noted  in  the  Introduction,  non-uniformity  in  fiber  spacing  introduces  a  new  aspect 
into  the  compressive  and  buckling  behavior  of  fiber  reinforced  composites,  namely 
transverse  internal  lateral  loads  associated  with  the  common  deformation  of  the  fibers. 
Following  the  statistics  of  spatially  distributed  data  and  the  concept  of  Voronoi  cell 
tessellation,  as  employed  to  represent  the  spatial  distribution  of  spherical  and  cylindrical 
inclusions  (Davy  and  Guild  (1988)),  we  assume  a  cumulative  distribution  function  for  the 
cell  size  2c  described  by  a  Poisson's  point  process 


P  (C  >  c)  =  exp  (-  2pc) 


10 


In  equati(xi.  (23)  ^  is  the  frequency  of  Voronoi  cells  in  a  unit  length,  with  a  mean  cell  size  of 
The  above  consideration  is  subject  to  the  restriction  that  fiber  regions  cannot  overitq), 
namely  oh  ("Gibbs  hard  core  process").  Therefore  equaticm  (23)  is  noodified  to  read 

P(C  >  c)  =  exp  (-  2^'  (c  -  h))  (24) 


Since  is  still  the  expected  value  of  the  Voronoi  ceU  size,  namely. 


one  obtains 


H-^=E(2c)  =  -|  2c^P(C>c)dc 

Ll  ss - — - 

^  l-2ph 


(25) 


Equations  (23)  -  (25)  can  be  expressed  in  terms  of  the  fiber  volume  fraction  Vf,  as  employed 
in  equation  (13).  Let  Vf  denote  the  average  ("nominal")  value  of  the  fiber  volume  fraction 
and  2c  =  |i"^  the  average  length  of  the  Voronoi  cells,  then  Vf  =  h/c  =  2hp.  Consequently, 
we  have  _ 

2h(l-Vf) 


and 


P  (C  >  c)  =  exp 


1-Vf'h  / 


Therefore,  the  cumulative  probability  that  the  fiber  volume  fraction  Vf  exceeds  a  value  Vf  is 


P(Vf>V,)=l-P(C>c)=l-enp[-^^(i-l) 

The  probability  density  distribution  which  ctxresponds  to  equation  (26)  is 

p(Vf)  =  --^p(Vf>Vf)  =  -^-Lcxpf — ^(^-ij 
Kv  ll  dVf  '  *  ^  1-VfV?  I  1-VflVf 


(26) 


(27) 
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ComputaticMud  results  for  p  (Vf)  vs.  Vf  are  shown  in  Figure  5  for  three  nominal  (average) 
values  of  Vf  (Vf »  0.3,  0.6,  and  0.9). 


The  Compressive  Response  with  Randomly  Spaced  Fibers 

The  probability  density  P  (^f)  given  in  equation  (27)  was  incorporated  into  the  formu¬ 
lation  expressed  in  equations  (11)  and  (13)  and  employed  to  predict  the  compressive 
response  of  Gr/PEEK  (APC-2)  composite  with  Vf  »  0.6  at  a  temperature  of  T  »  2VC. 
Based  upon  the  data  of  Guynn  et  al.  (1992),  the  nonlinear  shear  behavior  of  the  PEEK  resin 
was  fitted  by  a  Ramberg-Osgood  expressicm 


iA  , 


(28) 


where  GS*  ~  1.3  GPa  as  in  the  previous  section,  A  =  94.4  MPa  and  n  =  0.12.  In  addition, 
we  took  e  =  bJL  =  1/200  as  before  and  assumed,  somewhat  arbitrarily,  resin  failure  to  occur 
at  YilV  =  Yu  ~  10%.  The  latter  assumption  was  guided  by  the  observed  tensile  failure  at 
Eu  -  4%-5%  for  PEEK  at  room  temperature  reported  by  Johnston  et  al.  (1991).  The  shear 
stress-strain  response  considered  in  the  foregoing  representation  is  shown  in  Bgure  6. 

The  solution  to  equation  (13),  with  Y(0)  -  ^(2)  together  with  (27)  and  (28) 

was  obtained  numerically.  Note  that  equation  (28)  was  supplemented  by  xjPy  =:  0  for  lygy  I  > 
Yu  •  To  implement  the  numerical  solution,  the  field  equation  (13)  was  expressed  by  finite 
differences  as  given  by  Na  (1979),  and  solved  it^tively  by  a  quasi-linearization  method. 

In  the  above  implementation,  the  probability  distribution  function  of  the  Voronoi 
cells,  ^Vf),  was  evaluated  at  100  equally  spaced,  discrete  values  of  Vf  varying  between  Vf  - 
0  and  Vf  s  1.0.  With  the  exception  of  Hgures  12  and  13,  all  computations  were  performt  1 

for  Vf  =  0.6. 

Further  details  of  the  numerical  schemes  are  given  in  the  Appendix. 

Upon  attaining  convergence  to  a  prescribed  degree  of  accuracy,  the  computational 
program  gives  the  values  of  v^  Y,  Y'  arul  Y',  as  well  as  the  shortening  of  the  coluixu  A. 
Results  for  the  non-dimensionalized  lateral  deflection  vVL  and  for  the  slope  Y  vs.  X  are 
shown  in  Figures  7  and  8  for  three  values  of  non-dimensional  compressive  loads  X,  namely 
X  =  10, 20  and  26.4.  The  latter  value  corresponds  to  the  buckling  load,  since  no  equilibrium 
configuration  could  be  computed  for  X>26.4.  The  variation  of  ySV>  strain  in  the 

matrix,  vs.  the  distance  X  at  X;=26.4  is  shown  by  the  solid  line  in  Figure  9.  This  variadon  is 
contrasted  with  the  variation  of  ySV  X  for  uniformly  spaced  fibers  at  the  same  load  level. 
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as  shown  by  the  dashed  line,  and  against  the  variation  of  ySV  vs.  X  for  uniformly  spaced 
fibers  at  X=29.S,  which  is  the  maximal  load  level  attained  in  the  uniformly  spaced  case,  as 
shown  by  the  dotted  Une.  All  the  plots  in  Hgure  9  correspond  to  Vf  =  0.6  (in  the  case  of 
random  spacing  Vf «  0.6  and  the  results  are  plotted  for  the  cell  widi  Vf  =  0.6). 

Substitution  of  the  numerically  obtained  solution  for  v^  into  equation  (6)  determines 
the  lateral  load  q(x)  for  each  Voronoi  cell,  as  specified  by  its  fiber  volume  fraction  Vf. 
Results  for  q  vs.  the  non-dimensional  distance  X  =  x/L  are  shown  in  Figure  10  for  a  typical 
"matrix  rich"  cell,  with  Vf »  0.25,  at  load  levels  corresponding  to  X  »  10, 20  and  the  buck¬ 
ling  value  X  s  26.4.  Similar  plots  are  shown  in  Figure  1 1  for  a  "matrix  poor"  Voronoi  cell, 
with  Vf  s  0.95.  Note  that  sufficiently  low  levels  of  X,  i.e.  X  =  10,  yield  small  values  of  lat¬ 
eral  load  q,  while  increasing  levels  of  X  raise  the  magnitude  of  q.  It  is  especially  interesting 
to  note  the  "spikes"  in  the  plots  of  q  vs.  X.  These  localized  amplifications  occur  at  places 
where  ySV  attains  its  ultimate  value  at  some  Voronoi  cell,  with  the  sharpest  spike  located 
near  the  place  where  lysy  I  =  Yu  at  the  Voronoi  cell  under  consideration.  For  instance,  the 
spikes  in  q(X)  for  X  =  20  in  Figure  10  occur  at  X  =  0.15  and  X  =  0.3,  which  are  the  loca¬ 
tions  where  lysy  I  =  Yu  at  the  Voronoi  ceUs  of  fiber  volume  fractions  Vf  =  0.99  and  Vf  = 
0.98,  respectively,  at  X  =  20.  (Obviously,  the  matrix  material  in  those  cells  failed  over  the 
ranges  of  0.15  <  X  <  0.5  and  0.3  <  X  <0.5,  respectively).  On  the  other  hand,  the  sharp 
spike  at  X  s  0.25  for  X  s  26.4  in  Figure  1 1  is  associated  with  ySV  attaining  its  ultimate  value 
Yu  within  the  very  same  Voronoi  cell  (with  Vf  =  0.95)  considered  in  that  figure,  while  the 
remaining  peaks  are  associated  with  shear  failures  in  other  cells.  Peaks  which  occur  at  loca¬ 
tions  X  <  0.25  are  due  to  failures  in  cells  with  values  of  Vf  >  0.95,  while  spikes  located  at  X 
>  0.25  are  due  to  failures  within  more  resin-rich  Voronoi  cells.'^ 

Comparison  between  Figures  10  and  1 1  shows  that  resin-rich  Voronoi  cells  are  sub¬ 
jected  to  relatively  lower  lateral  loads.  This  observation  is  attributable  to  the  fact  that  the 
above  mentioned  cells  sustain  shear  strains  YSy  of  comparatively  smaller  magnimdes. 

Predicted  axial-stress  axial-strain  relations  and  compressive  strengths  under  mono- 
tonically  increasing  compressive  loads  are  illustrated  in  Figure  12  for  various  values  of  Vf. 
The  continuous  lines,  terminating  at  points  which  corresponds  to  failure,  correspond  to 
uniformly  spaced  fibers,  while  symbols  represent  computational  results  for  the  case  of 

,1/2 

tit  may  seem  that  lateral  equilibrium  is  not  satisfied  for  the  individual  Voronoi  ceils  since  j  q(X)dX^0  in 

the  plots  shown  in  Figures  10  and  11.  However,  due  to  symmetry  about  X-0  and  X  =  05,  J  q(X)dX 
indeed  vanishes. 
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randomly  spaced  Voronoi  cells  with  filled  symbols  representing  failure.  The  stress-strain 
responses  shown  in  Figure  12  are  dtxninated  by  the  last  term  on  die  tight  side  of  equation  (4) 
and  thus  remain  nearly  linear  until  failure.  In  Figure  13,  predicted  levels  of  cmnpressive 
strength  are  plotted  versus  fiber  volume  firaction,  Vf,  for  uniformly  and  randomly  spaced 
fibers.  Note  that  random  spacings  yield  lower  values  of  compressive  strength  and  suggest  a 
linear  relation  between  strength  and  Vf,  which  accords  with  experimentally  observed  trends 
by  Piggott  and  Harris  (1980). 

Figures  14  (a,  b)  exhibit  plots  of  fiber  curvatures  versus  non-dimensional  distance  X 
at  various  levels  of  non-dimensional  compressive  load  X.  Note  the  significant  increase  in 
curvature  for  the  randomly  spaced  case  (Figure  14  (b)),  as  compared  with  the  uniformly 
spaced  case  (Figure  14  (a)).  If,  according  to  Yin  (1992),  Idnks  occur  when  fibers'  curvature 
attains  a  critical  value,  then  Figures  14  suggest  that  random  spacing  yield  Idnks  at  lower  load 
levels. 

Unlike  the  circumstance  of  uniformly  spaced  fibers  with  bi-linear  shear  response  of 
the  matrix,  the  computational  scheme  for  randomly  spaced  fibers  carmot  be  extended  to  pre¬ 
dict  post-buckling  behavior  such  as  shown  in  Figure  2.  The  specific  values  of  the  computed 
compressive  failure  stresses  are  listed  in  Table  1.  That  Table  exhibits  the  effects  of  the 
nominal  volume  fraction  Vf,  the  amplitude  of  geometric  imperfection  5o/L,  and  the  presence 
or  absence  of  an  ultimate  value  of  matrix  shear  strain  Yu-  The  Table  also  illuminates  the  effea 
of  random  fiber  spacing. 

3.  CONCLUDING  REMARKS 

This  article  presented  a  mechanics  model  for  the  compressive  response  and  failure  of 
uni-directioiudly  reinforced  polymeric  composites  loaded  parallel  to  the  fiber  direction.  The 
model  accounted  for  the  non-linear  shear  response  of  the  resin,  including  its  ultimate  shear 
strain,  and  incorporated  two  kinds  of  geometric  imperfections,  namely,  initial  fiber  waviness 
and  random  fiber  spacings.  Heretofore,  the  latter  kind  of  imperfection  has  not  been 
considered  elsewhere. 

The  non-linear  response  of  the  matrix  was  accounted  for  by  means  of  the  non-linear 
field  equation  (6)  for  the  lateral  displacement  v^.  In  general,  the  above  equation  could  be 
solved  numerically  up  to  failure.  Nevertheless,  in  some  special  circumstances,  it  was 
possible  to  generate  a  solution  into  the  post  buckling  range. 

Both  kinds  of  geometric  imperfections,  initial  fiber  waviness  and  random  fiber  spac¬ 
ings,  were  shown  to  substantially  reduce  the  compressive  strength  of  the  composite.  How¬ 
ever,  random  fiber  spacings,  when  combined  with  the  foregoing  non-linear  shear  response 
of  the  matrix,  was  shown  to  introduce  imbalances  in  the  support  furnished  by  the  matrix 
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against  fiber  microbuckling  —  resulting  in  highly  localized  internal  transverse  loads  on  the 
fibers.  The  emergence  of  these  transverse  loads  alludes  to  the  possibility  of  transition  from 
microbuckling  to  microkinking  of  the  deformed  fibers.  However,  it  is  impossible  to  explme 
this  matter  any  further  within  the  context  of  the  BemouUi-Euler  beam  theory  utilized  in  the 
present  article  since  this  theory  cannot  account  for  discontinuous  shear  deformations  within 
the  fibers.  Such  discontinuities  are  likely  to  occur  at  locations  where  the  matrix  reaches  its 
ultimate  strength  and  ceases  to  support  the  fibers,  and  the  highly  concentrated  transverse 
loads  predicted  by  the  present  analysis  reflect  the  indetermiruu:y  inherent  in  the  Bemoulli- 
Euler  theory  in  addressing  shear  response. 

A  remedy  to  the  above  inadequacy  may  be  found  by  emplo'  shear-deformation 
models,  such  as  the  Timoshenko  beam  theory,  to  represent  the  resp  )f  the  fibers.  This 
approach  was  employed  recently  by  Qiung  and  Weitsman  (1993),  where  it  was  shown  that 
random  fiber  spacing  indeed  causes  discontinuities  in  the  shear  strains  within  the  fibers 
These  discontinuities  indicate  the  emergence  of  kinks. 
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The  nonlinear  second  order  differential  equation  (13)  can  be  expressed  as 

Y"=<3(X,Y)  (al) 

where  the  prime  denotes  derivatives  with  respect  to  X,  and 

Q(X.Y)  =  I  p(Vf)  aHVf)  (1-Vri  F(j^ )  dVf-  >?Y  -  X^Yo 

An  enor  quantity  at  i-th  iteration  step  is  defined  as 

<|>^*)  =  Y"<‘>-Q(X,Y<‘)) 

Consequently,  upon  employing  a  Taylor  series  expansion,  the  subsequent  error  quantity  is 
given  by 

-  Y<‘>)  +  -  Y’*<‘>)  (a2) 

Notingthatj^l  =  - ^  and  =1,  we  obtain,  upon  imposing  =  0  in 

equation  (a2) 

Y"(i+i)  _  [^j^‘ =  CKX,Y<*))  -  (a3) 

Expression  (a3)  is  a  linear  ordinary  differential  equation  for  Y<*'^i)  involving  the  known 
results  of  the  previous  iteration  Y^>).  Note  that  the  derivative  of  Q  with  respect  to  Y  is 

^=1  p(V,)oWt)  f(^)  dV,-X2 

Furthermore,  upon  employment  of  the  Ramberg-Osgood  model,  we  have 


F  = 


X 


,  G?  /  „ 


Obviously,  the  boundary  conditions  in  equation  (15)  must  be  satisfied  in  every  iteration  step. 

The  linear  differential  equation  (a3)  is  solved  by  a  finite  difference  scheme  as 
follows.  Divide  the  abscissa  0<X<1  into  N  equal  intervals  of  length  h=l/N.  Then  at  each 
node  XsXn=nh  the  second  derivative  Y"  is  expressed  as 
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Using  the  above  relation,  equation  (a3)  can  be  converted  to  an  algebraic  equation  of 
the  form 

(a4) 

Here, 

jQ(X,„Y<‘>)-(^]fi’Y<») 

The  boundary  conditions  in  finite  difference  scheme  are  =  0  and 

The  system  of  equations  (a4)  can  be  represented  as 

y^(i+i)Y(i+i)  =  (a5) 

where 


"bf'> 

1 

0 

1 

1 

1  b?"'^ 
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1 
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Y^i+t) 

\’‘N 

4-‘>l 

Equation  (a5)  can  be  solved  by  means  of  the  LU  decomposition  (Na,  1979).  Accordingly, 

the  matrix  A0+*)  is  decomposed  into  the  product 

Here, 


L<i+i)  = 


-p(M, 

0 

1 

0 

1  Pf" 

1  y(‘^'> 

U<i+u  = 
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and 


p(W) . 

(n=U...N-l) 

p(W) .  _  Y«*i)  („  .  2,3...N-1) 

pS^'>.C’-2yS:1’ 


Denoting 

2(i*i)  _  jjO+i)  y(»+»)  (a6) 

equation  (a5)  is  transfcnmed  to  where  the  components  of  are 

computed  by 

=  („  .  2.3...N-1) 

4r‘'=(4»‘>-2z<jr.‘>)/pg^» 

The  recursive  relations  between  and  are  obtained  from  equation  (a6)  as 

vCt+i)  _  ,(i+i) 

‘N  “*N 

=  (n=  N-lJJ-2 1) 

The  values  of  express  the  solution  to  equation  (13)  at  the  (i+l)th  iterative  step.  When 


attains  a  constant  value  within  a  prescribed  tolerance,  the  iteration  is  halted 

n^l 

and  post-processed  to  compute  deflection,  shear  strain  and  stress,  lateral  stress  and  other 
quantities. 
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Figure  1. 


(a)  A  Hber  composite  modelled  as  a  two  dimensional  lamellar  region 
consisting  of  fiber  and  matrix  plates;  (b)  a  defonned  single  cell 


Figure  2.  The  scaled  compressive  displacement  A/L  at  X=O.S  vs.  applied  compressive 
stress  Cc  for  various  fiber  voluix»  fractions  Vf  (symbol  conesponds  to 
the  circumstance  of 

Figure  3.  The  dimensionless  length,  ^  =  0  J  of  the  inelastic  zone  of  matrix  shear 

response  ([ySVl  ^  Yy  in  equation  (18))  vs.  iq)plied  compressive  stress  for  Vf 
=  0.3, 0.6  and  0.9. 

Figure  4.  The  variation  of  the  matrix  shear  strain  ySV  non-dimensionalized 

distance  X  along  the  fiber/matrix  interface  at  several  values  of  non- 

dimensionalized  applied  compressive  stress  X.  Fiber  volume  fraction  Vf  » 
0.6.  Onset  of  departure  from  linear  elastic  matrix  shear  response  at  X  Xy  = 
30.79,  maximum  compressive  stress  at  X  =  Xmu  =  30.81. 

Figure  5.  Distribution  of  local  fiber  volume  fraction  for  randomly  spaced  fiber 
conqtosites  with  average  fiber  volume  fraction,  Vf,  of  0.3,  0.6  and  0.9. 

Figure  6.  Shear  constitutive  relation  of  PEEK  at  21*0  based  on  Guynn's  (1992) 
estimation  with  shear  failure  strain  assumed  at  10%. 

Figure  7.  Non-dimensionalized  deflection,  vVL,  vs.  X  for  randomly  spaced  fiber 
conqxjsite  with  Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10, 

20  and  26.4.  Failure  shear  strain  Yu  is  10%,  and  X  =  26.4  is  the  compressive 
strength  of  the  composite. 

Figure  8.  Solution  Y  of  the  governing  equation  for  randomly  spaced  fiber  composite 
with  Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10,  20  and 

26.4.  Failure  shear  strain  Yu  is  10%,  and  X  =  26.4  is  the  compressive 
strength  of  the  composite. 

Figure  9.  Comparison  between  the  matrix  shear  strain  within  the  Voronoi  cell  with  Vf = 
0.6  in  randomly  spaced  fiber  composite  under  its  failure  load  X  =  26.4  and 
the  matrix  shear  strain  for  uniformly  spaced  fiber  composite  under  the  same 

load  level  as  well  as  with  its  own  failure  load  X  =  29.5.  Vf  is  0.6  for  both 
cases  (RS  and  US  designate  randomly  and  uniformly  spaced  fiber  composite, 
respectively). 
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Figure  10. 

Figure  11. 

Figure  12. 

Figure  13. 
Figure  14. 


Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf=  0.25  in  randomly 
spaced  fiber  composite  with  Vf  -  0.6  at  various  levels  of  non-dimensional 

compressive  loads  X.  The  load  X  »  26.4  corresponds  to  the  failure  strength 
of  the  conqrosite. 

Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.95  in  randorrrly 
spaced  fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional 
compressive  loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength 
of  the  composite. 

Dimensionless  displacement  -A/L  at  X=0.5  vs.  applied  compressive  stress. 
(Solid  lines  are  for  uniformly  spaced  fiber  composite.  Symbols  are  for 
randomly  spaced  fiber  composite.  The  ends  of  lines  and  the  filled  symbols 
indicate  compressive  failure  strength  for  uniform  and  random  spacings, 
respecdvely.) 

Compressive  strengths  of  uniformly  and  randomly  spaced  fiber  composites 
vs.  fiber  volume  fiction. 

Pre-buckling  Curvature  of  fiber  layer  in  the  case  of  (a)  Uniformly  spaced 

fiber  composite  and  (b)  Randomly  spaced  fiber  composite.  In  both  cases  Vf 

=  0.6. 
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(a) 

_  q(x) 


(b) 

Hgure  1.  (a)  A  fiber  reinforced  composite  modelled  as  a  two  dimensional 
lamrilar  region  consisting  of  fiber  and  matrix  plates;  (b)  a  deformed  single  celL 
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Figure  2.  The  scaled  compressive  displacement  d/L  at  X=0.5  vs. 
applied  compressive  stress  Oc  for  various  fiber  volume  fiactions  Vf 

(symbol "+"  corresponds  to  the  circumstance  of  |yg‘y(X=^)j  =  Yy). 


Figure  3.  The  dimensionless  length,  ^  =  O.S  -  of  the  inelastic  zone 
of  matrix  shear  response  ([Y!V|  >  Yy  in  equation  (18))  vs. 
iqrplied  compressive  stress  for  Vf  =  0.3, 0.6  and  0.9. 
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Figure  4.  The  variation  of  the  matrix  shear  strain  ySV  vs.  the  non-dimensionalized 
distance  X  along  the  fiber/mairix  interface  at  several  values  of 

non-dimensionalized  applied  coii^)ressive  stress  X.  Fiber  volume  fraction  Vf  =  0.6. 
Onset  of  departure  from  linear  elastic  matrix  shear  response  at  X  =  Xy  =  30.79, 
maximum  compressive  stress  at  X  =  X^u  ==  30.8 1. 


Figure  5.  Distribution  of  local  fiber  volume  fraction  for  randomly  spaced 
fiber  composites  with  average  fiber  volume  fraction,  Vf,  of  0.3, 0.6  and  0.9. 
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Hgure6.  Shear  constitutive  relation  of  PEEK  at  2rC  based  on 
Guynn's  (1992)  estimation  with  shear  failure  strain  assumed  at  10%. 


Hgure  7.  Non-dimensionalized  deflection,  vVL,  vs.  X  for  randomly 
spaced  fiber  composite  with  Vf  =  0.6,  under  compressive  loads 
corresponding  to  X  s  10, 20  and  26.4.  Failure  shear  strain  Yu  is  10%, 
and  X  =  26.4  is  the  compressive  strength  of  the  composite. 
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Figure  8.  Solution  Y  of  the  governing  equation  for  randomly  spaced  fiber  composite 
with  Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10, 20  and  26.4. 
Failure  shear  strain  Yu  is  10%,  and  X  =  26.4  is  the  con^ressive  strength  of  the  composite. 


Vf  =  0.6  in  randomly  spaced  fiber  composite  under  its  failure  load  X  =  26.4  and 
the  matrix  shear  strain  for  uniformly  spaced  fiber  composite  under  the  same  load 

level  as  well  as  with  its  own  failure  load  X  =  29.S.  Vf  is  0.6  for  both  cases 
(RS  and  US  designate  randomly  and  uniformly  spaced  fibCT  composite,  respectively). 
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Figure  10.  Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.25  in  randomly 
spaced  fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensitmal  conqnessive 
loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength  of  the  composite. 


Figure  11.  Lateral  stress  q(X)  vs.  X  on  a  Vcaondi  cell  with  Vf = 0.9S  in  randomly 
spaced  fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional  compressive 
loads  X.  The  load  X  »  26.4  corresponds  to  the  failure  strength  of  the  composite. 
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Figure  12.  Dimensionless  displacement  -A/L  at  X=O.S  vs.  applied  compressive 
stress.  (Solid  lines  are  for  uniformly  spaced  fiber  composite.  Symbols  are  for 
randomly  spaced  fiber  con^site.  The  ends  of  lines  and  the  filled  symbols 
indicate  conq)ressive  failure  strength  for  uniform  and  random  spacings,  respectively.) 


Figure  13.  (Compressive  strengths  of  uniformly  and 
randomly  spaced  fiber  composites  vs.  fiber  volume  fraction. 


29 


Hgure  14.  Pire-buckling  Curvature  of  fiber  layer  in  the  case  of  (a)  Uniformly  spaced 
fiber  composite  and  (b)  Randomly  spaced  fiber  composite.  In  both  cases  Vf  «  0.6. 
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